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Melting in Na„ clusters described with an empirical embedded-atom potential has been reexam- 
ined in the size range 55 < n < 147 with a special attention at sizes close to 130. Contrary to 
previous findings, premelting efi^ects are also present at such medium sizes, and they turn out to 
be even stronger than the melting process itself for Naiaa or Naias. These results indicate that the 
empirical potential is qualitatively unadequate to model sodium clusters. 



Introduction 

The early days of cluster thermodynamics were first 
mostly concerned with theoretical of numerical studies 
of simple rare-gas aggregates. Beyond the now generally 
accepted idea that melting in finite atomic clusters ap- 
pears as a first-order transition rounded by size effects fi 
it was also observed in simulations that these system can 
exhibit dynamical coexistence,^ a process in which the 
cluster fluctuates in time between its solidlike and liq- 
uidlike states. 

While no experiment has yet permitted to test these 
predictions on the very same clusters, recent observations 
achieved in the group lead by Haberland^*^ have pro- 
vided valuable qualitative and quantitative information 
about the way sodium clusters melt, by extracting the full 
caloric curves from photodissociation measurements. In 
particular these authors reported that the melting point 
and the latent heat of fusion both vary strongly and non- 
montonically with sizcs^ These results significantly at- 
tracted the attention of theoreticians, who since then 
tried to interpret or reproduce these unexpected com- 
plex variations using a variety of models&^iSiSiiSiiiiiS 
Further evidence was also seeked in indicators, which 
are alternative to the caloric curves, such as the elec- 
tric polarizability,-'^- the photoabsorption spectrum,^** or 
the ionization potential. 

Up to now, none of the above cited theoretical works 
has been able to reach a fully satisfactory quantitative 
agreement with experiments. The situation is largely 
due to the expected strong interplay between geomet- 
ric and electronic effects, which could be responsible for 
the variations of the thermodynamic quantities. How- 
ever, the previous simulations have brought some clues 
about the relevance of the various models and potentials 
used to describe simple metal clusters. For instance, it 
was seen that the distance dependent tight-binding (TB) 
model developed by Poteau and Spiegelmani^ overesti- 
mated melting points by more than 20 percents. This 
was interpreted as the consequence of the parameteriza- 
tion of this model, carried out on small clusters only but 
not on bulk propertiesii On the other hand, both the em- 
pirical embedded-atom model (EAM) potcntialSiii^iS 
and orbital-free density functional^ calculations lead to 
a notable underestimation. Unfortunately, more realistic 
simulations still lack the extent of phase space sampling 



required for a precise computation of the melting point 
above 50 atoms. 



Numerical experiments 

In our previous works, we concluded that melting in 
sodium clusters actually occurs differently in the small- 
est and in medium to large sizes clustersi^ While the 
caloric curve of the smallest clusters (having less than 
about 80 atoms) exhibit several features due to multi- 
ple isomerizations prior melting, the heat capacity of the 
larger sizes mainly has one main peak indicating melt- 
ing in a non ambiguous way. This could partly explain 
why experimental measurements do not find such pre- 
melting features above 55 atoms, and also why they have 
difficulties in getting a clear, single-peak picture below 
this size. More recently we noticed that premelting ef- 
fects could also be artefactually due to poorly converged 
simulationsi^ 

Based on our previous results, the disagreement be- 
tween experimental measurements and the results of both 
the TB and EAM models seemed mostly quantitative. 
One could hope in getting a much better agreement by 
suitably modifying the parameters, after including both 
molecular and bulk properties, possibly through allowing 
these parameters to become size-dependent. However, 
the range of sizes that was investigated by us and by 
others was quite limited^SiLiiiiS with only very few sizes 
above 80. Moreover, the mediocre agreement for the la- 
tent heats of fusion lead us revisit the problem with newly 
available simulation methods and less ambiguous tools of 
analysis. 

We have performed exchange Monte Carlo canonical 
simulationsiSi of the clusters Na„ with 55 < n < 145, 
n being a multiple of 5 in this range, plus the following 
sizes n = 59, 93, 127, 129-131, 133, 139, 142, and 147. 
The clusters are again described using the same empirical 
many-body EAM potential whose parameters are given 
in Ref. T?, and each simulation consisted of 10^ cycles fol- 
lowing 2 X 10^ equilibration cycles for each of the 31 tra- 
jectories characterized by their temperature Ti = i x 10 K 
for 1 < i < 30 plus i = 1.5. The starting structures for 
all clusters was always chosen to be the result of basin- 
hopping global optimization carried out with 10 sets of 
10^ quenches. All were found to be based on the icosahe- 
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dral motif. The absence of any different structural mo- 
tif such as octahedral or decahedral for the sizes studied 
here is very favorable for the simulations to reach equilib- 
rium and not fall into broken ergodicity problems, as one 
major cause for such difficulties precisely lies in the en- 
ergy landscape having several funnels»i2*i^ Each cluster 
was simulated 5 times independently, with different ran- 
dom seeds, and we used a hard-wall spherical container 
with radius -Rmax = 7n^/^. The caloric curves were con- 
structed from the distributions of potential energies using 
a multihistogram technique. In Ref. |3, we calculated the 
latent heat of melting, L, as the integral of the heat ca- 
pacity minus the Dulong-Petit contribution. This lead 
to appreciable overestimates, mostly due to the neglect 
of anharmonicities. Here we proceed similarly to the ex- 
perimental approach of Schmidt and coworkers,"* namely 
by fitting the low- and high-temperature parts of the in- 
ternal energy as straight lines and defining L as the gap 
between these lines at the melting point Tmcit- The melt- 
ing point itself is defined as the temperature at which the 
last heat capacity peak has its maximum: in cases where 

[k) 

there are several peaks each centered around T^dt 
true melting temperature is taken as in.ax.k{Tl^^i^ \ . The 
low- and high-temperature parts are defined as T < 50 K 
and T > 250 K, respectively. Therefore, if premelting 
events are present between these limits, they will con- 
tribute to the latent heat. 



Results and discussion 
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FIG. 1: Melting point Tmcit and latent heat of melting L of 
sodium cluster clusters versus their size. The open squares 
are the experimental results of Schmidt et alfi-, the full circles 
are for Monte Carlo simulations using the empirical potential. 

The variations of Tmcit and L with size are sketched in 
Fig.^along with the latest experimental data of Schmidt 



and coworkers.^ Except for very low sizes, we first no- 
tice that the melting points computed here are signifi- 
cantly lower than the ones we previously reported^ This 
is an unfortunate consequence of using the q-jumping 
method with un inappropriate Tsallis parameter. How- 
ever they are comparable to the ones obtained by Garzon 
and coworkers in the microcanonical ensemble. Since 
we expect the differences between this ensemble and the 
canonical ensemble to become smaller and smaller for 
increasing sizes, this agreement appear as a mutual con- 
firmation of the convergence of both calculations, even 
though they rely on very different numerical experiments. 

Interestingly, the experimental and theoretical varia- 
tions of Tineit with size appear to be related to each other, 
especially if we shift the simulation data around size 120 
by about 15 atoms and 70 K. This had not been noticed 
before, and is probably fortuitous. But it may also hide 
that some mechanisms causing the strong variations in 
the experimental results close to 140 atoms are indeed 
the same here, only occuring sooner. 

In Fig. Hone should also notice that the melting point 
at n = 55 is not the highest, as clusters having 59 or 60 
atoms are more resistant to an increase in temperature. 
However, we were not able to extract any latent heat for 
these sizes (nor for n = 65 and 70) due to very broad heat 
capacity peaks. Therefore one should maybe not give too 
much importance to the melting points extracted from 
these curves. 

While the melting temperatures are usually well be- 
low the experimental data (except in the vicinity of 130 
atoms), the computed latent heats show a reasonable 
overall agreement, even though the complex variations 
are not as sharply seen as in the measurements of Schmidt 
et al.^ A definitive comparison would require one to ex- 
tend the range of sizes. As far as latent heats are con- 
cerned, the difference between the present results and our 
previously published datai comes nearly entirely from the 
different way of estimating L, which is now much closer 
to the experimental way. 

The heat capacity curves in the range 125 < n < 135 
are all plotted in Fig. |2 versus temperature. Despite 
strong changes from one size to another, a regular evo- 
lution can be seen from 125 atoms and above this size. 
The heat capacity consists of two peaks, the melting (or 
high-temperature) peak being centered near 203 K. The 
smaller peak, denoted as premelting peak in the follow- 
ing, goes from 100 to about 180 K in a quite continuous 
fashion. The premelting peak is surprisingly strong, and 
can be clearly distinguished from the melting peak. Strik- 
ingly, it is even stronger than the melting peak itself at 
size 133. 

To interpret these curves, we have chosen to focus on 
Nai33, by carrying periodic quenches from instantaneous 
configurations extracted from the Monte Carlo simula- 
tions for all trajectories. We thus gathered nearly 14000 
different isomers. Fig. |2| shows the energy of these iso- 
mers versus their rank, as well as the discrete spectrum 
of isomers versus the temperature of the trajectory from 
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FIG. 2: Heat capacities of sodium clusters calculated with 
exchange Monte Carlo using the empirical potential, in the 
size range 125 < n < 135. The inset shows the curve obtained 
for Nai33 using the tight-binding quantum model. 



which they were quenched. Both graphs clearly indicate 
a correlation between the repartition of isomers in energy 
space, their number and the heat capacities. Small vari- 
ations in Cy occur below 150 K, but are hardly visible 
when compared to the main peaks. They are related to 
the first few hundreds of isomers, which involve the mi- 
gration of a limited number of some of the missing atoms 
on the external layer of the icosahcdron (the 135-atom 
cluster has Ih symmetry with all vertex atoms missing). 

The presence of two major peaks in Cy is consistent 
with the two main increases in the number of isomers 
having less than a given energy. Most new isomers ap- 
pearing between 150 and 190 K have their rank between 
1000 and 5000 in the upper part of Fig. |31 Looking at 
their structure reveals that they are all still based on the 
two-layer icosahedron, but that the third layer is hardly 
recognisable. Hence this case of premelting is an extreme 
illustration of surface melting following preliminar sur- 
face reconstruction.--'^ Eventually, above 200 K the icosa- 
hedral structure is completely lost and the true (volume) 
melting takes place. 

If now we look more closely at the results obtained for 
Nai35 J we notice that the main peak is indeed located at 
the same temperature as the premelting peak in Nai33, 
and that the true melting peak of the latter cluster has 
been replaced by a right shoulder. Thus the same phe- 
nomena seem to be present in Naias, only with different 
relative magnitude. This case could be referred to as 
'post-melting'. 

To some extent, this progressive evolution of the pre- 
melting and melting peaks can be compared to what has 
been observed in small Lennard- Jones (LJ) clusters by 
Frantz.^^ In these systems, a premelting peak starts ap- 




5000 1 0000 

Isomer 




100 200 
Temperature (K) 



FIG. 3: Quenching analysis of the Monte Carlo trajectories for 
Nai33. Lower panel; spectra of isomers versus temperature. 
Upper panel: energy versus isomer rank. 



pearing close to the size 31 corresponding to the compe- 
tition between Mackay-type and anti-Mackay-type icosa- 
hedral structures. The premelting peak remains as the 
cluster size increases, and it is shifted to higher tempera- 
tures before becoming higher than the melting peak itself 
near the size 38 Because of the structural similarities 
between LJ clusters and the present sodium clusters, it 
is likely that the present observations express the same 
qualitative mechanisms. However, a significant quanti- 
tative difference can be seen in the caloric curves of LJ 
clusters and sodium clusters, as the melting peak is never 
really well resolved for van der Waals systems until it has 
replaced the premelting peak. Here the two peaks are 
of very similar widths, but their respective heights vary 
strongly, their total contribution to the latent heat being 
nearly constant. 

Up to now, experimental data did not find any evi- 
dence for any pronounced premelting peak in the heat 
capacities of charged sodium clusters. ^■'^'^ This indicates 
that the present empirical potential is qualitatively unad- 
equate to describe these systems. In particular, beyond 
a simple scaling of the parameters of the potential, we 
do not expect the use of explicit, size-dependent param- 
eters to improve the situation notably. We repeated the 
above simulation for Naiaa using the more realistic quan- 
tum distance-dependent tight-binding (TB) Hamiltonian 
described in Ref. 0, but we had to reduce the statistics 
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to 10^ cycles following 2 x 10^ equilibration cycles (per 
trajectory) for the computation to be tractable. Even 
though melting points were shown to be overestimated fi 
we also noticed^ that premelting effects were quite re- 
duced using this model, in better consistency with ex- 
periments. The heat capacity computed from exchange 
Monte Carlo simulations is reported in the inset of Fig.|21 
For this calculation we neglected the (weak) effects of 
nonzero electronic temperaturei^ The starting configura- 
tion was taken as the same one as for the classical poten- 
tial, but we did not find any more stable structure during 
the course of the simulation. 

At first sight, the caloric curve looks similar to the 
classical result, with a clear premelting peak. However 
using the TB Hamiltonian has two consequences. First, 
the premelting peak is much lower than the melting peak, 
which is in agreement with our previous general observa- 
tion that the empirical potential overemphasizes premelt- 
ing features. More importantly, premelting also occurs 
much closer in temperature to the melting peak itself. 
This also suggest that premelting is not seen experimen- 
tally simply because it is too broad. 

Conclusion 

In the present work, we have obtained some evidence 
that premelting effects in the caloric curves of sodium 
clusters could be present at unexpectingly large sizes. 
We also found that some clusters could exhibit 'post- 
melting', a process in which the premelting effect is 
stronger than the actual melting peak. In the cases stud- 



ied here, these effects seem to be associated with sur- 
face reconstruction of the third icosahedral layer, and 
thus seem to be of character similar to what occurs in 
Lennard-Jones clusters having about 35 atomSfSS, 

One consequence of the above results is that explicit 
empirical potentials are not fully reliable for predicting 
melting points in small sodium clusters, not only because 
they do not allow one to reproduce the complex variations 
observed by the group of Haberland,* but mostly because 
they exhibit prominent premelting peaks not seen in ex- 
periments. 

Calculations performed using the quantum tight- 
binding model also predict a premelting phenomenon 
near 133 atoms, but the corresponding anomaly of the 
heat capacity is much smaller than the melting peak, as 
well as closer to it. In this respect, it resembles more the 
measurements by Schmidt et alA 

Even though our calculations overestimate premelting 
effects, they provide insight into the possible causes for 
the nonmonotonic variations of the melting point. In 
particular, they suggest that such variations may reflect 
premelting becoming actual melting. The discrepancies 
with the present work would then be ascribed to a pos- 
sible merging of the premelting feature into one shoulder 
of the melting peak, but not necessarily on the low tem- 
perature side. 
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